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Abstract 

The Steady State (SS) probability distribution is an important quantity needed to characterize the steady state 
behavior of many stochastic biochemical networks. In this paper, we propose an efficient and accurate approach to 
calculating an approximate SS probability distribution from solution of the Chemical Master Equation (CME) under 
the assumption of the existence of a unique deterministic SS of the system. To find the approximate solution to 
the CME, a truncated state-space representation is used to reduce the state-space of the system and translate it to 
a finite dimension. The subsequent ill-posed eigenvalue problem of a linear system for the finite state-space can 
be converted to a well-posed system of linear equations and solved. The proposed strategy yields efficient and 
accurate estimation of noise in stochastic biochemical systems. To demonstrate the approach, we applied the 
method to characterize the noise behavior of a set of biochemical networks of ligand-receptor interactions for 
Bone Morphogenetic Protein (BMP) signaling. We found that recruitment of type II receptors during the receptor 
oligomerization by itself doesn't not tend to lower noise in receptor signaling, but regulation by a secreted co- 
factor may provide a substantial improvement in signaling relative to noise. The steady state probability 
approximation method shortened the time necessary to calculate the probability distributions compared to earlier 
approaches, such as Gillespie's Stochastic Simulation Algorithm (SSA) while maintaining high accuracy. 



Introduction 

Many biological networks exhibit stochasticity due to a 
combinatorial effect of low molecular concentrations 
and slow system dynamics. One important biological 
context where stochastic events likely have a large 
impact is the Bone Morphogenetic Protein (BMP) sig- 
naling pathway. BMPs make up the largest subfamily of 
the Transforming Growth Factor-/? superfamily and are 
involved in numerous processes including growth, dif- 
ferentiation and diseases [1]. Due to their potency at 
driving development, they are also of great value for 
stem-cell differentiations in cell culture. BMPs activate 
near maximal signaling at InM concentration, have very 
slow binding kinetics and require oligomerization 



* Correspondence: dumulisiapurdue.edu 

'Department of Agricultural and Biological Engineering, Purdue University, 
West Lafayette, USA 

Full list of author information is available at the end of the article 



between multiple receptor subunits [1]. These properties 
naturally lead to conditions for significant and long- 
duration stochastic fluctuations in cellular signaling. 
Interestingly, variability of BMP signaling appears to be 
very low in vivo, while it is very high in stem cell culture 
studies [2]. To understand the differences between 
in vivo and in vitro signaling and determine how various 
receptor oligomerization events might alter the signal 
and noise, a more efficient means of solving the steady 
state distributions for stochastic model was needed that 
would allow for continuation of both parameters and 
levels of the BMP pathway components. 

Stochastic regulation can negatively impact the robust- 
ness of the system [3,4] or instead, constructively con- 
tribute to the phenotypic variation [5-7] in a species. 
In stochastic reaction networks, the state of a species tra- 
verses different trajectories in a probabilistic manner and 
the distributions of states can be difficult to predict. 
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As more biological data is available, stochastic modeling is 
becoming increasingly popular to estimate properties in 
networks where the time evolution of the system is unpre- 
dictable and dependent on unavoidable randomness inher- 
ent to the system. The complete solution can be calculated 
from a Chemical Master Equation (CME) [8-10], that is 
based on a Markovian approach that captures the inherent 
randomness of biochemical systems. 

The Chemical Master Equation (CME) describes the 
dynamics of the probability distribution of a species of 
chemical reactions. Precisely, the CME captures the rate 
of change of probability that a system will be in state X at 
time t for all the species of the system. Solution of the 
CME is practically intractable due to the curse of dimen- 
sionality, as the state-space of the system becomes enor- 
mously large with increases in the species number and 
concentrations (number of states n N , for N — > species, 
n — > copies of each species). Moreover, the system often 
involves interactions between different time-scales (slow 
and fast reactions, frequent and infrequent transitions 
between states) [11], which add further complexity. 
Instead, numerical approaches are commonly used 
[12-14] to realize the CME of the stochastic system. 

In the analysis of stochastic biochemical networks, 
steady state probability distributions for each species in 
the system are determined to measure variability about 
the deterministic steady state. The deviation around the 
solution contributes to stochastic noise that can be quan- 
tified by measuring the coefficient of variation A = ^ 
(the ratio between the standard deviation a and the mean 
level of species concentration y) obtained by solving the 
CME [9]. 

Frequently, Monte-carlo based simulation approaches 
[9,13] are used to solve stochastic problems. But, there 
are drawbacks in this approach for networks where the 
dynamics of different intermediate states of the system 
are unknown and continuation of several parameters is 
necessary to explore the system's dependency. As a 
screen of parameter values becomes necessary in such a 
scenario, the Monte-carlo based approach doesn't prove 
to be efficient, as it generally takes longer time to 
numerically simulate the process and satisfy the imposed 
conditions. Moreover, simulation times increase with 
increases in the total number of molecules, species and 
the number of interactions between species. 

In order to ameliorate the computational cost and 
complexity, we present a method here to approximate 
the steady state probability distribution by 1) reducing 
the system's state-space to a finite dimension using 
truncated state-space method [15] and 2) subsequently, 
translating an eigenvalue problem associated with a 
CME to a system of linear equations. We illustrate that 
the eigenvector (for an eigenvalue = 0) that represents 



the steady state probability distribution can be solved 
algebraically by approximating it as a system of linear 
equations. Previously, the influence of stochastic fluctua- 
tions on system behavior was studied also in [16], where 
a moment closure approach was applied to reduce the 
complexity associated with the identification of distribu- 
tions. In contrast to the previous studies, here we use a 
truncated state-space for steady state probability distri- 
bution approximation, which is arguably more general 
since we make no assumptions on the relationships of 
the moments of the distribution. 

The usefulness of the proposed method is illustrated 
by considering the example networks of BMP signaling, 
as described earlier in [17]. Here we examine two poten- 
tial mechanisms of BMP signaling: 1) regulation between 
the type I and type II receptors, and 2) regulation by 
secreted co-factors, such as Crossveinless-2 (Cv-2). First, 
we apply the approach to the recruitment of Type II 
receptor into a BMP-bound type I receptor complex to 
see if such step of receptor oligomerization contributes 
qualitatively to the noise profile of the system. Following 
this, we extend our earlier work that focused on extra- 
cellular regulation of BMP signaling by SBPs to evaluate 
the calculation approach and compare results to the 
Type I/Type II receptor system. 

Unlike SBPs, which tend to improve the signal to 
noise ratio, we did not see a significant change in sto- 
chastic variability with inclusion of Type II receptors. 
This result supports a previous assumption made in [4] 
where the recruitment of Type II receptor was excluded 
to simplify the modeling steps while characterizing the 
noise profile of a SBP regulated BMP signaling system. 

Methods 

Proposed approximation method 

The Chemical Master Equation (CME), which is a set of 
first order differential (ODE) equations, demonstrates 
loss and gain of probabilities of discrete states of a system 
[10] and is often applicable to represent the stochasticity 
of the system. Consider a well-stirred system at thermal 
equilibrium of N different species {Si, S 2 , ■ • ■ , S N } with 
{X u X 2 , ■ ■ ■ , X N } molecules respectively, participating in 
a total of M biochemical reactions R M , where fi - 1, 2, . . . 
M. The state of such a system is represented by the copy 
number (X n molecules) of each species (S„) at any given 
time t and is represented as X = [Xi(t), X 2 (t), . . . , Xjjf)]. 
Unless a non-zero initial state is assigned, the default 
initial species concentrations are always zero (X„(t = 0) = 
0, where 1 < n < N). 

Two other quantities are further required to con- 
struct the system: 1) a state-change vector and 2) 
propensity functions [8,12,14,18] for the reactions R^, 
fi. = {1, 2, . . . , M}. The state-change vector for 
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reactions is defined as = [v lfi , 
where v nfi represents the change in 
species S n , caused by the occurrence 
the underlying biochemical system, 
fully define the system and the time 
probability function P(X, t) can be 
solution of the Chemical Master 
[8,14,18]: 



v 2fi> • • • , V NfJ \ , 

concentration of 
of Rf, reaction of 
These equations 
evolution of the 
obtained by the 
Equation(CME) 



(1) 



a„(X)P(X,t)) [(X + iv) e ^] 



Here, [(X +v fi ) e O] is 1 if X +v fl e Cl and 0 other- 
wise. The CME representing the rate of change of prob- 
ability P(X, t) in an in finitely large state-space X e O is 
given by taking Cl to be the non-truncated space: Cl = 
N N , N = {0, 1, 2 . . .} 

In Eq.l, represents the propensity function to 
account for transition from a given state X to any other 
state, and indicates the stoichiometry of the reaction 

that results in such a transition. Eq.l is a linear system 
of differential equations and may be rewritten as follows: 



d_P 
eft 



LP 



(2) 



where P is the probability distribution P(X, t) for a 
vector X = [X lt X 2 , . . . , X N ] and L is the time indepen- 
dent connection operator. For the steady state (SS) dis- 
tribution P ss , we have: 



i)Pss X >0; VXefi ii)^P ss (X) = l 



Xe£2 



(3) 



and in) LP SS = 0, 



We assume that the deterministic steady state (SS) is 
unique. The non-truncated state-space Cl can be 
replaced with a truncated state-space [15,19] to 
approximate the probability distribution P(X, t). We 
define the truncated space as: 



£2 = {X : «j < Xi ; < fii, Vi} 



(4) 



where <*,- and /3, are extendable left and right boundaries 
of the truncated state-space. This approach is similar to 
that in [20], in which it is shown that the approximation 
based on the truncated space converges to the true steady 
state distribution as the limits of the truncated state-space 
converge to the limits of the original space. 

The truncated state-space representation implies that 
given some e > 0, for a sufficiently large /3, > 0 and suffi- 
ciently small a, > 0, the steady state probability distribu- 
tion P ss (X) is approximated to within e: 



£> S (X) = 1-s 



For an optimal SS probability approximation, e should 
be made as small as possible. In the truncated state- 
space, Eq.3(iii) is represented as: 



LP« = 0 



(5) 



where £ is a matrix of propensities in q . To get 
the entries in £ we use Eq.l modified so that 
P (X, £) = 0 if X <=£ Q anda M (X) = 0 if X + v tJi £ Q . In the 
truncated state-space Q , Eq.5 is an eigenvalue problem for 
eigenvalue X = 0 and the eigenvector p ss can be obtained 
algebraically, or with an iterative algorithm for a large, 
sparse matrix £ . 

Instead of finding the eigenvector, which can be an ill- 
conditioned problem when there are nonzero eigenvalues 
close to 0, we translate the problem to a well-conditioned 
system of linear equations as follows. 

We first evaluate the deterministic steady state (Y 0 ) of 
the system, and then select state X 0 of the discrete system 
closest to this deterministic steady state, where X 0 = 
round{Y Q ). Taking p ss to be the solution of Eq. 5 and 
using the fact that the deterministic steady state solution 
is unique, we observe that P ss (X 0 ) is among the largest 
entries of p ss . The states in are labeled as 1, 2, . . . , K 
with state X 0 denoted by /. 

Then 



ss ml 



P 

1 ss 



[pI...P ss ...p? s ] t /Hs 



(6) 



[<7i,...,4,-_i, l,^ + i ...cj K ] T 



pk 

1 SS 



where q k = k = 1 . . . K, cj,= \. With 

Pis 

cj = [cf\, 1, ...%] Eq.5 now becomes Lq = 0- Let Li, be 
the ^ column of l . Expanding Lcj by column and rear- 
ranging gives the following well-conditioned problem: 



L k q k = -Lj or, 

L'cj' = -Lj 



(7) 



In Eq.7, £' is the matrix £ with column / removed and 
c( is cj with entry cjj removed. The error criterion for 
the system is checked for the calculation of p ss until a 
satisfactory value is obtained (see algorithm 1 for further 
details). 

Application 

In order to demonstrate the usability of the proposed 
steady state probability approximation method, we present 
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here two example networks (example network 1 and 2) 
from Bone Morphogenetic Protein(BMP) mediated signal- 
ing, and characterize the stochastic behavior of the sys- 
tems. In the example network 1, we consider the role of a 
specific extracellular protein, Crossveineless-2(Cv-2), 
which is part of a class of proteins known as Surface-asso- 
ciated BMP-binding Proteins (SBPs) [4]. Cv-2 has the abil- 
ity to regulate the stochastic noise in BMP signaling, and 
in this example we demonstrate that the role of Cv-2 is 
heavily dependent on reaction kinetics of the network: for 
some sets of parameter values, Cv-2 increases the coeffi- 
cient of variation of the steady state signaling distribution, 
while for other parameter values it decreases the coeffi- 
cient of variation. 

In the second example network, we consider a model 
simplification strategy as used in [4,17]. This strategy is 
to omit a Type II receptor recruitment step from the 
receptor oligomerization in a BMP patterning model, 
under the assumption that the simplification step does 
not affect the outcome of a BMP-mediated patterning 
model. The obtained results by the use of steady state 
probability approximation method provide a numerical 
justification for the aforementioned simplification. 

Background 

During embryonic development, positional information is 
transduced by morphogens to underlying cells that 
respond to the concentration gradient of morphogen and 
eventually differentiate into distinct cell types [21]. For 
example, Decapentaplegic (Dpp), a drosophila homologue 
of BMP2/4, forms a spatial profile to pattern dorsal tissues 
in Drosophila development [21], In a canonical BMP sig- 
naling pathway, dimeric ligands bind to receptors and 
form a heterotetrameric complex that consists of two 
Type I and two Type II receptors. The heterotetrameric 
receptor complex then phosphorylates the intracellular 
signal transducer Smad and the phosphorylated Smad 
forms a complex with a co-Smad. Subsequently, the 
Smad/Co-Smad complex accumulates in the nucleus and 
regulates gene expressions in a concentration dependent 
manner [17,22]. 

BMP regulation occurs at many points along the path- 
way, and a lot of focus has been on identifying and 
understanding how the ligand activity is regulated in the 
extracellular environment by secreted binding proteins. 
These include molecules such as Cv-2, HSPGs, among 
other reviewed in [1]. A focus of this work is to gain a 
better understanding of how regulation in the extracellu- 
lar region impacts cell signaling noise, and eventually 
cell-to-cell variability. 

Example networks 

In many biochemical networks, where dynamics of the 
intermediate interactions of different species (proteins) 



and molecular complexes are largely unknown, screen- 
ing plays a significant role in the classification of 
dynamics-dependent network behavior. For example: 

1. In a biochemical network where a class of secreted, 
surface-associated BMP binding proteins (SBPs) such 
as, Crossveinless-2 (Cv-2, node D as in Figure 1) [4] 
is allowed to regulate BMP signaling, the intermediate 
dynamics of the system that result in the formation 
and decoupling of a transient state BMP:Type I:Cv-2 
(node M as in Figure 1) are largely unknown. 

2. In the patterning modeling of BMP signaling path- 
ways, it is often argued as a simplification strategy that 
omitting the step of recruitment of a Type II receptor 
to a bound BMP:Type I receptor complex doesn't 
affect the outcomes of patterning models [4,17,23,24]. 
While valid in the deterministic sense, it is not clear 
how this reduction impacts our estimates for noise in 
the sytem. 

In these systems, we apply our SS probability approxi- 
mation method to evaluate the probability distribution of 
different species and calculate the mean (^<), standard 
deviation (cr) and the coefficient of variation (A = ^ ; 
defined as the ratio between the standard deviation and 
the mean of any species) of the species distribution. 
Together with this information, we can screen the net- 
work for largely unknown dynamics of the intermediate 
interactions and classify solutions according to a model's 
ability to meet specific performance objectives. 

BMP-signaling regulation by SBPs 
Signaling network 

The single-cell local stochastic model that includes extra- 
cellular BMP(A), receptors (B), and SBPs such as, Cv-2 
(D) with biochemical interactions, rate parameters, and 
connectivity is based on the network shown in Figure 1. 
Mass balance equations are listed below: 

Rt-.tfhA R 2 : A 0 R 3 :A + Bhc 
R 4 : C A + B R 5 :A + D% E 
R 6 : E -H A + D R 7 :B + E%Z 
R 8 : Z -H B + E R g :C + D^Z 
R 10 : Z -H C + D 

Out of all complexes (C = BMP:Type I Receptor = BR, E 
= BMP:Cv-2, Z = BMP:Type I receptor: Cv-2), available 
experimental evidence suggests that only ligand-bound 
receptors C (BMP:Type I Receptor = BR) initiate signaling 
to regulate downstream gene expression [17]. To focus on 
the noise compensation by regulation of receptors by 
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Figure 1 Example network of BMP signaling. BMP signaling is mediated by BMPType I Receptor (C) forms either by direct interaction 
between BMP (A) and Type I (B) or via an intermediate state with BMPType l:Cv-2 (Z). Initially B interacts with A (Type I receptor), D (Cv-2) and 
forms C (BMPType I receptor) and E (BMP:Cv-2) complexes respectively. Then the complexes C, E can generate the intermediate state BMP:Cv-2: 
Type I complex (node Z) by interacting with D and B respectively. In the network, only BR(node C) has the ability to turn on downstream 
signaling. Upon BR formation, the complex recruits type II receptor, and later initiates the phosphorylation of intracellular Smad protein. 
Signaling leads to pSmad accumulation within the nucleus and gene expressions of BMP targets. 



SBPs, the extracellular level of BMP (A) is treated as a 
parameter A = ^ and the interactions (1-10) are simpli- 
fied accordingly. For example, in the simplified model 
reaction R 3 is represented R' • B — > C • 

The simplified model as obtained from reactions Ri to 
R 10 , has 5 species {Si, S 2 , ■ ■ ■ , S 5 } and is described com- 
pletely by a total of 8 different chemical reactions. Time 
evolution of all species quantities are specified by a state 
vector X(t) = [X^t), X 2 (t), . . . , X 5 (t)] T and state-change 
vector (fi - 1, 2, . . . 8), corresponding to all reactions 
that describe the system. For example, when fi - 1, V\ is 
[-1 +1 0 0 0] T for reaction R\ of the simplified system. 
In this example network, techniques like those in [25] 
were used to verify that there is a unique steady state 
equilibrium and this ensures the applicability of the 
algorithm for this example network. Numerically, we 
used the polynomial root finding package hom4ps2 to 
ensure that there was only one equilibrium in the posi- 
tive orthant [26]. It's worthwhile to mention that a simi- 
lar approach is adopted in example network 2 to ensure 
the unique deterministic steady state. For both the net- 
works, we numerically determined the deterministic 
steady state value Y 0 using Newton's method as 



incorporated in SBTOOLBOX2 [27]. A generalized algo- 
rithm for simulation according to the steady state 
approximation as outlined in Methods section is given 
in algorithm 1. 

Algorithm 1 Evaluate steady state (SS) distribution: 
LP SS = 0 

Require: Unique deterministic SS solution X 0 

1. Reaction Networks with N Reaction R lt . . . , R N 

2. Choose: e, y 0 , y 

3. Solve: ODE for steady state(SS) = Y 0 and find dis- 
crete state X 0 closest to Y 0 , where X 0 = round(Y 0 ). 

4. Initiate, a h /?,•; where a; = (X 0 )j - y 0 , 
Pi = {X 0 )i + yo 

5. Determine: Truncated state-space Q as shown in 
Eq.4 and £ as described after Eq.5 

6. Determine: Column j of £ corresponding to X 0 

7. Form fj and Lj as describe after Eq. 7. 

8. Solve: L'cj' = -Lj 



Find 



where 



cj' =[cj 1 ,...,q j - 1 ,cj j+1 ,...cl K ] T and f) > 0 and 
r\ = 1 + <\\ is chosen so that ^ SS ( X ) = 1 

1=1,1)0 X€€l 
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10. Verify: 

if E p ss (X) > e for 8 . = a or §. = p. then 

Xef2,X, = (5, 

at m ai-Y 

Return to 5 
end if 

In the algorithm, the values of Jo, J are problem 
dependent based on the anticipated spread of the steady 
state (SS) distribution. Larger y and y 0 favor a larger q , 
with correspondingly better accuracy, but this comes at 
the expense of a larger state-space and more time 
required to solve the Eq.7. The parameter s also controls 
the accuracy of the solution. In the truncated state- 
space, the tail of the distribution is essentially pushed in 
to the main part of the distribution and smaller e means 
that less of the tail is changed. 
Simulation and discussion 

The binding kinetics between BMPs (species A, Figure 1) + 
receptors (species B), and BMPs + Cv-2 (species C) are 
largely known from the biacore analysis data [28,29]. How- 
ever, the kinetic data associated with the intermediate tri- 
partite complex BMP:Cv-2:Type I receptor (species Z, 
Figure 1) are currently unknown. In order to better under- 
stand the dependence of the steady state distribution on 
the kinetic parameters, we performed a parameter screen 
for the forward and reverse reaction rates (k ±s , s = 3, 4) for 
the formation and decoupling of species Z. For each of 
these four parameters, we use 5 evenly-spaced points on a 
logarithmic scale with the range [10 1 to 10 1 ] wAf V 1 for 
the forward rates and [10" to 10 ] 5" for the reverse reac- 
tion rates. This produces a parameter grid that contains a 
total of 625 different parameter vectors. 

For an appropriate comparison of the noise attenua- 
tion both in the presence and in the absence of Cv - 2, 
species C (BMP: Type I Receptor = BR) concentration 
should remain the same regardless of the intermediate 



dynamics. During simulation, the amount of available 
receptors (B) was fixed at 100, and a maximum of 30% 
receptor occupancy was allowed for the screening of the 
network. To ensure an equal amount of BR formation 
for each parameter vector, we modified the level of free 
ligand (A). For computational tractability, the screen is 
limited to a maximum continuation of 200 Cv-2 mole- 
cules, which allowed us to capture responses for all 625 
different parameter sets. 

In order to quantify noise in the system we measured 
the coefficient of variation (A = — ) that relates the stan- 
dard deviation (a) to the mean (//) level of bound recep- 
tors. The parameter screen on the intermediate 
dynamics yields three primary qualitative subclasses for 
Cv-2 behavior in regulation of extracellular BR (C) fluc- 
tuation amplitude: i) reduced amplitude ii) increased 
amplitude and Hi) mixed amplitude behavior [4]. Three 
primary types of responses of Cv-2 action on BMP fluc- 
tuations are shown in Figure 2a-c. As seen from Figure 
2a, Cv-2 leads to a reduction of BR noise amplitude, 
and this is true for all Cv-2 e [0,200]. The value of the 
coefficient of variation (A) decreases for both increases 
in the level of bound receptor and the level of Cv-2. 
The subclass of increased amplitude demonstrates that 
increasing the level of Cv-2 in the system increases the 
value of the coefficient of variation (A Cv2 * o > A Cv2 = o) 
and is valid for the range of Cv-2 values considered in 
the screen (for a detailed discussion on this, interested 
readers can refer [4]). Lastly, mixed amplitude is classi- 
fied as type Hi, which demonstrates that Cv-2 can both 
increase and decrease the level of stochastic noise in the 
system (Figure 2c). 

To clarify Cv-2 action further, we calculated per- 
centage change in the amplitude using 

(Acii-2=105~Aq 



(% noise change = 



x 100) for each 



parameter set output for a given amount of BR produc- 
tion (30 complexes) and Cv-2 value (105 molecules). 




Cv-? , o 
02 = 0 



200 



100 



10 20 
Number ol BR 



10 20 

Numb* •• W 



10 20 

N Mb) ' ' ' BR 



30 



5,o 



~l~n-f-i 



Now* change (%) 

Figure 2 Screening result, a, b, c) Coefficient of variation (A = for BR (C) (1 to 30 molecules) formation is shown both in the presence of 
Cv - 2 = 105 molecules (gray) and in the absence of Cv-2 (black), where a, b, c represents the noise attenuation, noise amplification and biphasic 
subclasses of Cv-2 behavior respectively, d) Histogram of noise change for all 625 parameter sets demonstrates that the majority of solutions 
result in noise attenuation. 
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Based on the percent change of the coefficient of varia- 
tion, we classify the screening outcome: negative percent 
change corresponds to noise attenuation whereas the 
positive change gives noise amplification. A histogram 
of all 625 parameters are shown in Figure 2d and in [4]. 
The implication of the screening result is that Cv-2 
clearly reduces the variability of receptor activation 
throughout the range of Cv-2 tested. However, as 
demonstrated in Figure 2d, such a phenomenon as 
exhibited by SBPs like Cv-2 is found to be highly para- 
meter dependent. 

During the simulation, the kinetic rate constants for 
the intermediate complex BMP:receptor:Cv-2 (node Z, 
Figure 1) formation and decoupling were chosen from 
the parameter grid and representative kinetic rate con- 
stants for three different type of Cv-2 responses are 
enumerated in Table 1. 
Analysis of Type II receptor recruitment process 
In the signaling network shown in Figure 3, recruitment 
of Type II (= Ri) receptors can happen in two different 
ways: 1) BMP binds with Type I (= RJ first and subse- 
quently, recruits Type II receptors to form a tripartite 
complex BMP:Type I: Type II (BR 1 R 2 ), and 2) BMP 



directly interacts with Type I and Type II separately, 
and an intermediate state forms a tripartite BMP:Type I: 
Type II complex. In both situations, BMP:Type I:Type II 
complex (B7? 1 7? 2 ) is the sole signaling complex responsi- 
ble for the activation of downstream pathways. 

All possible biochemical interactions that represent 
the ligand binding with Type I receptors and further 
recruitment of Type II receptors are: 

ri : B + i?j % BRi r 2 : BKi B + i?i r 3 : B + R 2 % BR 2 

r i :BR 1 k ^>B + R 2 r 5 : Bi?i + R 2 ^ B«iR 2 r 6 : BRiR 2 -H Bfti + R 2 

r 7 : BR 2 + J? x BRiR 2 r 8 : Bi?ifi 2 -H BR 2 + i? r 

The chemical interaction of Case II can easily be 
obtained from the interactions {r^ to r 8 ) of Case III 
(Figure 3) by equaling the kinetic rate constants k ±2 and 
k+4 of Case III to zero. For the kinetics, we relied on the 
published data [1]. The rate at which a Type II receptor 
is recruited upon formation of a BMP:Type I complex 
(BRi) is comparatively faster than the rate of BMPs and 
Type I receptors interactions [17]. However, exact values 
of the rates of formulation of (BR 1 R 2 ) complex are not 
readily available, and hence, parameters were screened 



Table 1 Kinetic rates. Figure 2(a,b,c) 

Figure k 3 (molecule" 1 sec" 1 ) k_ 3 (sec 1 ) It, (molecule" 1 sec" 1 ) k_ 4 (sec 1 ) 

2a 1.3282 0.0100 1.3282 0.0100 

2b 0.0133 1.0000 0.1328 0.0100 

2c 0.1328 1.0000 0.4200 1.0000 



Case I 



BMP 



Y Type I 

Y Type II 



Case II 





r y 

9 + T=t *=Y + Y 

Figure 3 Network cases for Type II recruitment analysis in context of Drosophila melanogaster development. Case I) Recruitment of Type 
I is overlooked here and it imitates the simplified model used in previous studies. In this type of network, BMPType I complex (Bfi,) acts as the 
sole signaling complex. Case II) Upon the formation of a BMP: Type I complex, subsequent recruitment of Type II receptor is considered here. But a 
direct interaction between BMP and Type II receptor doesn't happen in the network. Here, a tripartite complex BMPType IType 11 (Bft^) activates 
the downstream pathways. Case III) Similar to Case II, but with an exception that a direct interaction between BMP and Type II receptor is allowed 
to form a BMPType II complex (Bft 2 ) by changing Bft, to BR 2 . The kinetic equations are equivalent to the SBP system investigated in [4]. 
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over the physiological ranges with values between [10" to 
10 1 ] nM' 1 s' 1 for the forward rates and [10~ 3 to 10°] s 1 
for the reverse reaction rates. 
Simulation and discussion 

To simulate the networks (as shown in Figure 3) for the 
calculation of the coefficient of variation A, we applied 
the truncated state-space approximation. During the 
simulation, a target of 1 to 30 signaling complexes (BRi 
for Case I and BR±R 2 for Case II, Case III) in the extra- 
cellular region is considered so a direct comparison can 
be made for the coefficient of variation (A = j^) 
between B7?! and BR^. 



The coefficient of variation (A) for BR X R 2 remains very 
close to the coefficient of variation of BR 1 as shown in 
Figure 4a. Proximity in the coefficient of variation 
between BR 1 and BR l R 2 (as shown in Figure 4a) demon- 
strates that the stochastic variability of the system is not 
affected by the recruitment of the Type II receptor. It is 
also found that increasing the concentration of R 2 brings 
the coefficient of variation of BR^R 2 into very close 
agreement with the coefficient of variation of BR 1 as 
shown in Figure 4b. A similar outcome is obtained from 
the simulation of Case III of the Figure 3 and the result 
is shown in Figure 4c. Finally, all the outcomes are 



a 

It 



0 5 



BR1 
BR1R2 



No of molecules 




R2 = 50 
R2 ■ 100 
R2 = 1000 
BR1 



No of BR1R2 



SO 



C 
15 



15 



BR1 
BR1R2 





10 20 30 
No at molecules 



No ofBR1R2MoMciMt 



M 



Figure 4 Comparison of A. a) The coefficient of variation of Bftj (calculated from Case I Figure 3)and BR^R 2 complexes (calculated from Case I 
Figure 3) is compared. The variability of the system seems to be invariant in the presence of Type II. b) Concentration dependency of A as a 
function of R 2 . c) Same as plot "a", however, direct interaction of BMP and Type II is allowed as in Case III, Figure 3. It's clear that the stochasticity 
of the system does not change over the range of values tested, d) Summary of BR^R 2 formation and its impact on signaling noise. 



II 

< 




Direct method 
ET -28000 hrs 
ET -2800 hrs 
ET -280 hrs 



Processing Time - 1 .7 min 
Processing Time - 34.5 min 



No. 



of BR 1 R 2 



Reference 
1000 times slower 



t 
n 

< 



30 



UQ 



0 



No. of BR 1 R 2 



30 



Figure 5 Comparison between SSA and Direct SS method, a) In Gillespie's method larger 'End Time' (ET) is required (which translates into a 
higher processing cost and time) to ensure the accuracy of outcome. Three different ET: 280 hrs, 2800 hrs, 28000 hrs are shown, b) Effect of 
kinetics associated with 8ft, interacting with R 2 . The steps of interactions are clearly shown in Case II of Figure 3. 
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summarized in Figure 4d, where it is shown that regard- 
less of the different cases as shown in Figure 3 the coef- 
ficient of variation (A) of BR^R 2 is approximately equal 
to that of BR^ 

Additionally, it is also found from the simulated result 
that the rate at which the BMP:Type I recruits Type II 
receptor (Case II in Figure 3) also decides the effect of 
Type II recruitment process on the stochastic variability 
of the system. With a comparatively slower rate, the 
coefficient of variation tends to oscillate as observed in 
Figure 5b. When the recruitment rate is slower than the 
formation rate of BMP:Type I complex, free Type II 
receptors fail to get frequent access to BMPs via the 
BMP:Type LType II tripartite complex, and can cause 
the concentration of BMP:Type LType II complex to 
oscillate more than the case with a comparatively faster 
dynamics. Thus, mitigating noise is not a natural output 
of receptor oligomerization + transudction and instead, 
requires another co-factor such as Cv-2 [4]. 

Benchmarking of Direct SS approximation method 

Carrying out large-scale stochastic simulation can be 
time consuming but calculation of the approximate solu- 
tion via a truncated state-space can greatly improve the 
speed. In order to show the performance improvement in 
terms of computational cost and time of direct SS 
approximation in the analysis of stochastic biochemical 
networks, we benchmarked the method by comparing it 
with Gillespie's stochastic simulation algorithm (SSA) 
method [9] for numerical calculations of stochastic bio- 
chemical networks. In the benchmarking, the processing 
time taken by each method was calculated based on the 
steps in the blue box as mentioned in the flow chart dia- 
gram (Figure 6). The sample problem was calculated for 
both methods on the same hardware and software config- 
uration: Processor. Intel(R) Xeon (R) CPU E5405, 2.00 
GHz (quad-core), RAM: 16 GB, SBTOOLBOX2 [27] and 
Matlab R2010a with SiMBiology 3.0. 

The processing time and computed A values for a target 
BR 1 R 2 = 20 for Case II, Figure 3, is provided in Table 2 to 
show the accuracy and time gain that can be obtained if 
the proposed direct SS distribution approximation method 
is used. Gillespie's SSA takes longer to generate an output 
that contains enough information to calculate the distribu- 
tion as compared to the time taken by Direct SS approxi- 
mation method. The problem becomes severe when 
continuation of a multiple parameters are necessary to 
explore the system's parameter dependency as done pre- 
viously in [4]. 

In Table 2, the term 'End time (ET) in Gillespie's SSA' 
corresponds to the amount of time the system dynamics 
were allowed to evolve. The accuracy of the Gillespie's 
SSA approach depends on the 'End time in Gillespie's 
SSA'(directly contributes to the processing time) set in the 



r \ 

Read parameter sets 

from grid 



Read 


target 


complex {BR 1 R 2 ) 


X 


o 



4 

Determine ligand 
needed 




Initial condition: SS 
concentration 



; ; 

Proposed Method/ 
Gillespie's SSA 



Output 

Figure 6 Benchmarking of Direct SS approximation method. 

Benchmarking of Direct SS approximation method. 

^ 

model simulation, and is shown clearly in Figure 5a and 
Table 2. Very low propensities require long simulation 
times in Gillespie's SSA due to the infrequency of events. 
Accuracy of Gillespie's method for the sample example 
increases as the 'End time in Gillespie's SSA' is increased. 
This large simulation time in turn directly impacts the 
processing time, resulting in a large computational cost to 
achieve the desired accuracy (Table 2). 

Conclusions 

In this study, we illustrate an approach of determining 
the steady state probability distribution efficiently to 
carry on continuation in multiple variables within a 
large-scale parameter screen. The approach is demon- 
strated further with a couple of applications, where we 
investigated 1) the dynamic dependency of a class of 
proteins, known as SBPs, in the regulation of BMP sig- 
naling, and 2) the binding of Type II receptor in BMP 
signaling. The results suggest that the recruitment of a 
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Table 2 Benchmarking: Gillespie's SSA and Direct SS approximation for a target BR^R 2 = 20 



Method 


End Time (ET) in Gillespie's SSA(hrs) 


Processing Time (sec) 


A 


Direct SS approximation 


Not Applicable 


0.4 -0.6 


0.1707 


Gillespie's SSA 


28000 


90-95 


0.1705 




2800 


8-10 


0.1878 




1390 


4-5 


0.2254 



type II receptor in BMP signaling doesn't affect the sto- 
chasticity of the system over the range of concentration 
and parameters investigated. Direct calculation of the SS 
probability distribution can be successfully applied to 
systems with a unique deterministic SS solution, and 
future work will investigate similar approaches for other 
biochemical systems. 
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